Sei’s drum!
Michael Großbach
Draft 2025-01-10
library(brms)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(kfigr)
library(knitr)
library(patchwork)The data set was read in as is, the columns subj,
arm, cond, series and
no were converted to factors, the column
filename was deleted. A new column was added to the
dataset, patient, based on the subj prefix
(‘P’ or ‘S’), with a ‘0’ coding for healthy controls, and a ‘1’ for
patients.
read.csv("features.csv") %>%
mutate(filename = NULL,
patient = factor(if_else(grepl("^P",
subj),
1L,
0L)),
series = factor(series),
subj = factor(subj),
arm = factor(arm),
cond = factor(cond),
no = factor(no)) ->
drum_beats1 Exploration
1.1 Instruction Condition, and Arm Used
In a first attempt to get an overview, the descriptors’ density
estimates were plotted, broken down by arm and
cond.
dep_vars <- names(drum_beats)[which(!(names(drum_beats) %in% c("subj", "arm", "cond", "no", "series", "patient")))]
n_vars <- length(dep_vars)
plot_lst <- vector("list", length = n_vars)
plt_cnt <- 1
for (dv in dep_vars) {
if (plt_cnt == n_vars) {
# plt_pos <- c(2, 0.5)
plt_pos <- "right"
} else {
plt_pos <- "none"
}
if (plt_cnt %% 4 == 1) {
ylab <- "density"
} else {
ylab <- ""
}
plot_lst[[plt_cnt]] <- ggplot(drum_beats,
aes(x = .data[[dv]],
color = .data[["cond"]])) +
geom_density(na.rm = TRUE) +
scale_color_colorblind() +
ylab(ylab) +
geom_rug(alpha = 0.25) +
theme(legend.position = plt_pos) +
facet_wrap(~ arm)
plt_cnt <- plt_cnt + 1
}
print(wrap_plots(plot_lst,
ncol = 4) +
plot_annotation(
tag_levels = "A"))
Fig. 1. Descriptor densities along with raw data points (rug ticks on x axes), broken down by arm (dominant [D] vs non-dominant [ND]) and instruction condition (controlled = C, normal = N).
Several peculiarities are immediately apparent from the plots in Fig. 1:
- There are no substantial differences in the density distributions between normal (N) and controlled (C) strokes
attDur,LAT, andattFlathave very few unique data values (see rug ticks at bottom of plots)- Looking at the dominant arm alone
(and excluding the vars, the densities have a trend toward a larger variance in the controlled condition (in other words: a less pointed distribution) with additional or larger humps in the tails, respectivelyattDur,LAT, andattFlatfor the reason just mentioned)
The lack of any clear-cut differences in the raw data between conditions suggests that it will be hard to find any differences by modeling.
The limited number of unique values in the attack-related
measures (e.g. attDur only has 35 unique values in
a total of 1102 observations across all subjects, conditions, and sides;
that’s only 3 percent!) suggests that rounding errors propagated through
the calculations. This is probably due to very similar, i.e.
highly automated, and short attack times combined with the given
sampling frequency, resulting in few data points, which in course of the
calculations result in the observed phenomenon. But I’m just guessing
here, as I do not know anything about these descriptors and how they are
calculated or interpreted. Judging by the range of values of,
i.e., attDur (2.59, 9.8ms, across both sides and
conditions) and a sampling frequency of 48 kHz, this leaves us at 48 *
(9.8 - 2.59) = 346 points to choose beginning and end of the attack.
Given the supposed highly automatized motor program used to initiate the
stroke, along with the laws of physics at play here (no pun intended),
it is not very surprising to see very few unique values. My limited
knowledge—or rather, my ignorance of anything sound-related—aside, from
a statistical standpoint these measures do not seem suited to describe
any differences between the experimental conditions investigated
here.
Given the very few data points in the attack phase, any
descriptor derived from such a short period of time (.e.g..
attSPL) cannot be judged as being stable in the sense of
being reproducible. Hence I suggest to drop all attack-derived
descriptors.
# drum_beats %>%
# mutate(attDur = NULL,
# LAT = NULL,
# attSPL = NULL,
# attSC = NULL,
# attFlat = NULL,
# attSpecFlat = NULL) ->
# drum_beatsThe increase in variance for the conditions N < C does not come as a surprise as I assume normal also means highly trained and thus automatized, whereas controlled involves less automatization and more ‘individualness’ both within and between subjects.
1.1.1 Conclusions
Looking at the above investigated descriptors and Danielsen et al. (2015), it seems reasonable to
limit the modeling attempts to totDur, totSPL,
totSC, and TC.
But Sofia wrote on 2020-07-09: “Francesco and I have discussed a bit related to descriptors and we want to concentrate on the “transient” period (although the name probably will change). Spectral Centroid (transSC) should be one, and I suggest transFlat for the other. Francesco, does that sound reasonable?”
On 2020-07-22 both Sofia and Francesco agreed upon
transSC, transFlat, and
transCrest as the probably most important response
variables to look at.
1.1.2 Francesco’s comments
Update:
- The lack of data points for
attFlatis fixed by solving a bug in the feature extraction: now the downsampling ratio for the envelope extraction is reduced, and makes the MIRToolbox algorithm less sensitive to frames with an rms value of 0 (which returns an incorrect value of 0, since we have a geometric mean at the numerator). Now we have a larger spread which allows high values (1 = peaky envelope, 0 = smooth/flat envelope). - The crest factor is now calculated for the separate phases.
Regarding the attack phase: I agree with the remarks regarding duration. Given the fact that we are not comparing different instruments, I wouldn’t have expected a large variation. This is not surprising if we consider that our system is changing only slightly (same rototom, same drumstick, same action, a bit different tuning across subjects): in fact, the perceived timbral differences are so small that we are in trouble guessing on the descriptors.
The sampling frequency is even lower (\(f_s = 44100\) Hz). Even if we had a higher sample rate which could reveal some discrepancies in the attack durations, we would have to prove that they are perceptually relevant.
Therefore, I am happy to discard attDur and
LAT (which is obviously a log-transformed duration, only
there for the sake of consistency with the literature).
I am a bit more in doubt when it comes to discarding all the attack descriptors. Even if what Michael says is true from a statistical point of view, we should still be able to discriminate timbre on short time windows due to the high temporal resolution of our hearing. Attack phase descriptors (with the same definition of attack that we are using, which is most likely not coincident with perceptual attack) are employed in Câmara et al. (2020), and the Oslo group has a paper under construction which analyzes drum sounds in a similar manner (see OsloPlots.png in the OneDrive folder).
I am worried that merely taking the overall descriptors into account would introduce a lot of unnecessary and perceptually catching information — mostly the tonal part of the signal, i.e. the drum ringing in the last part of the decay. That’s why Sofia and I are suggesting to look at what we could call “main energy” or “early decay” phase (i.e. from max peak to temporal centroid).
Would it be feasible to set up 4 different models (i.e. one for each phase), at least in the univariate version?
As for the descriptors to include: although TC is
employed in Danielsen et al. (2015), the
PDFs are even more similar. I would go for Dur (except
attack?), SPL, SC, and one between
Flat, SpecFlat or Crest.
My (informal and biased) listening tests tell me that, at least for some subjects, I hear a pattern going towards a harsher (controlled) vs smoother (normal) timbre exactly at the hit point, plus slightly less (controlled) or more pitch/amplitude fluctuation. Hopefully this could be catched by spectral centroid, specrtral/temporal flatness, or crest factor. This should be independent of SPL unless the subject misinterpreted the instructions, therefore SPL acts as a sort of control variable in our model.
1.1.3 Model Considerations
There are several points that need to be considered before making a decision regarding the type of modeling to be done in this study, (1) the sample size places restrictions on the external validity; (2) data from small samples can be better modeled when regularization is in place to ‘tame’ the estimates; (3) the hierarchical structure of the data (subjects played several trials with either their dominant or non-dominant arm under two conditions) suggests a multilevel analysis of the data which would, in addition to the Bayesian regularization via priors, also results in shrinkage of the estimates; (4) given the small distribution differences between the two experimental conditions, along with the
- Given the extremely small sample size, any attempt to fit the data using traditional statistics, a.k.a. null-hypothesis significance testing, would make us vulnerable to all kinds of criticism. I therefore suggest to use Bayesian statistics instead, as it allows us to estimate probability distributions of parameters rather than confidence intervals around point estimates, and thus embraces uncertainty in estimates.
- Additionally, Bayesian regression uses prior probability distributions to arrive at sensible estimates. These priors regularize estimates, or draw them toward zero, a desirable effect which has long been recognized even in traditional statistics (e.g. ridge regression and lasso; Tibshirani (1996)).
- ANOVA (or regression with grouping variables) requires each subject’s trials to be averaged, say, within conditions or groups (or both), to be able to assess treatment effects; this results in the individual variation within a subject to get lost while it might have been very informative to include it in the analysis. Multilevel (or hierarchical) modeling allows to include the entire data structure (Fig. 2) so that no information gets lost through averaging, but all variation (both individual and treatment-driven) propagates through the analysis and ends up in the final results, allowing for more realistic credibility margins around estimates.
- Variance in traditional statistics is considered to be fixed (think homoscedasticity in ANOVA, or the \(\epsilon\) in regression formulas such as \(y_i = \beta_0 + \beta_1x_i + \epsilon\)), whereas it is an estimated quantity in Bayesian statistics and therefore can vary between conditions (or groups or subjects etc.) and thus capture varying variance in different conditions, or groups etc.
# chunk intentionally empty Group: Patient Healthy Control
/ \ / \
Instruction: Normal Controlled Normal Controlled
| | \ / | \
Player: 1 2 ... n 1 2 ... n
/ \
Side dominant non-dominant
/ | \ / | \
Trial 1 .. p 1 .. pFig. 2. The multilevel structure of a data set should be reflected in the analysis.
1.1.4 Conclusions
Will start with a simple univariate model, add predictors and interactions, then a bivariate model, and finally a quadruple-variate model and see where this leads us.
1.2 Individual Responses
The humps in the density distributions in section Instruction Condition, and
Arm Used made me curious where they might originate from. So in the
following graph the density estimates of one of the descriptors,
transSC, are plotted broken down by subj. and
then also by cond, separately for patients and healthy
subjects.
ggplot(drum_beats, aes(transSC,
color = subj,
linetype = cond)) +
geom_density(alpha = 0.1) +
scale_color_colorblind() 
Fig. 3. transSC
density distribution of participants, broken down by instruction
condition.
It is obvious that some participants do not differ substantially between the normal and the controlled condition, whereas others do, and even markedly so. Additionally, there seems to be quite a spread of the centers of distributions across a wide range of values, suggesting very individual drum sounds.
The wide distribution of centers of mass between individuals made me want to further break down the plot.
ggplot(drum_beats, aes(transSC,
color = subj,
linetype = cond)) +
geom_density(alpha = 0.1) +
scale_color_colorblind() +
facet_wrap(~ patient)
Fig. 4. transSC
density distribution of participants, broken down by instruction
condition and group (patients and healthy subjects).
The above plot is interesting in that it seems to show that the four
healthy subjects were more uniform than the patients in their
transSC distributions and also, with the exception of S2,
had very similar transSC distributions for the normal and
the controlled conditions. In the patients, two had very similar
distributions in both conditions (P1 and P5), whereas the two others
showed differing results for the two conditions. So with regard to
Sofia’s hypothesis (ch. Modeling) I’d argue, at
least for transSC alone, having drummers play normal and
controlled strokes would not allow subjects to tell the difference in a
listening test. But maybe it would qualify as a screening test for
movement disorders in drummers. Just a random thought.
While breaking it down it occurred to me that looking at individual variation (that is, between series) might also be enlightening.
ggplot(drum_beats, aes(transSC,
color = cond,
group = series)) +
geom_density(alpha = 0.1) +
scale_color_colorblind() +
facet_wrap(~ subj, nrow = 2)
Fig. 5. Density plots of individual series.
And it was! The plot in Fig.
5 made visible that P1 and P5 had consistently very low transSC
values. P3 was consistent within the normal instruction condition, but
had, on average, higher values, and with a lot more variation, in the
controlled condition. P4 had more variation in both conditions, and
higher transSC values in both; in other words, P4 was consistantly bad.
(But then again: what do I know what bad is wrt
transSC!)
The healthy subjects showed comparable variation and centers of mass within conditions, and all but one (S2) also across conditions.
Comparing the two rows of panels in Fig. 5 reveals that healthy subjects have more variation than the patients.
Although tempting, we probably shouldn’t get carried away and generalize to the populations of healthy drummers and ones with movement disorders, respectively.
For completeness’ sake the same exploratory analysis was carried out
for transFlat and transCrest`:
ggplot(drum_beats, aes(transFlat,
color = subj,
linetype = cond)) +
geom_density(alpha = 0.1) +
scale_color_colorblind() 
Fig. 6.
transFlat density distribution of participants, broken
down by instruction condition.
ggplot(drum_beats, aes(transFlat,
color = subj,
linetype = cond)) +
geom_density(alpha = 0.1) +
scale_color_colorblind() +
facet_wrap(~ patient)
Fig. 7.
transFlat density distribution of participants, broken
down by instruction condition and group (patients and healthy
subjects).
ggplot(drum_beats, aes(transFlat,
color = cond,
group = series)) +
geom_density(alpha = 0.1) +
scale_color_colorblind() +
facet_wrap(~ subj, nrow = 2)
Fig. 8.
transFlat density plots of individual series.
ggplot(drum_beats, aes(transCrest,
color = subj,
linetype = cond)) +
geom_density(alpha = 0.1) +
scale_color_colorblind() 
Fig. 9.
transCrest density distribution of participants, broken
down by instruction condition.
ggplot(drum_beats, aes(transCrest,
color = subj,
linetype = cond)) +
geom_density(alpha = 0.1) +
scale_color_colorblind() +
facet_wrap(~ patient)
Fig. 10.
transCrest density distribution of participants, broken
down by instruction condition and group (patients and healthy
subjects).
ggplot(drum_beats, aes(transCrest,
color = cond,
group = series)) +
geom_density(alpha = 0.1) +
scale_color_colorblind() +
facet_wrap(~ subj, nrow = 2)
Fig. 11.
transCrest density plots of individual series.
1.3 Covariance Matrix
new_names <- drum_beats
names(new_names) <- gsub(pattern = "trans",
replacement = "earlyDecay",
x = names(new_names),
fixed = TRUE)
cov_mat <- cov(new_names[-c(1:5, ncol(new_names))],
use = "pairwise.complete.obs")
kable(cov_mat, digits = 2)| totDur | attDur | decDur | earlyDecayDur | LAT | totSPL | attSPL | decSPL | earlyDecaySPL | totSC | attSC | decSC | earlyDecaySC | TC | totFlat | attFlat | decFlat | earlyDecayFlat | totSpecFlat | attSpecFlat | decSpecFlat | earlyDecaySpecFlat | totCrest | attCrest | decCrest | earlyDecayCrest | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| totDur | 8580.54 | 32.18 | 7641.34 | 907.01 | 5.93 | -125.00 | -205.08 | -41.76 | -90.45 | -449.09 | -7391.47 | 1448.04 | -6574.23 | 999.35 | -0.79 | 1.26 | -2.41 | 2.57 | 1.05 | -0.17 | 1.31 | 0.12 | -90.78 | -13.10 | 27.04 | -39.39 |
| attDur | 32.18 | 1.15 | 30.22 | 0.81 | 0.23 | 0.26 | -0.66 | 0.25 | 0.43 | 8.36 | -171.34 | 47.75 | -101.81 | 2.68 | -0.02 | 0.02 | -0.02 | 0.01 | 0.01 | 0.00 | 0.01 | 0.00 | -1.26 | -0.26 | -0.03 | -0.40 |
| decDur | 7641.34 | 30.22 | 6838.06 | 773.06 | 5.66 | -108.21 | -177.56 | -36.46 | -76.60 | -325.99 | -7202.60 | 1404.24 | -5920.35 | 867.57 | -0.87 | 1.19 | -2.28 | 2.27 | 0.94 | -0.17 | 1.16 | 0.11 | -82.96 | -12.27 | 25.25 | -35.30 |
| earlyDecayDur | 907.01 | 0.81 | 773.06 | 133.13 | 0.04 | -17.05 | -26.86 | -5.55 | -14.28 | -131.46 | -17.53 | -3.96 | -552.06 | 129.09 | 0.09 | 0.05 | -0.10 | 0.28 | 0.11 | 0.01 | 0.13 | 0.02 | -6.57 | -0.57 | 1.83 | -3.68 |
| LAT | 5.93 | 0.23 | 5.66 | 0.04 | 0.05 | 0.10 | -0.08 | 0.10 | 0.14 | 1.11 | -37.72 | 8.53 | -19.80 | 0.40 | 0.00 | 0.00 | -0.01 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -0.26 | -0.05 | -0.01 | -0.08 |
| totSPL | -125.00 | 0.26 | -108.21 | -17.05 | 0.10 | 35.41 | 35.05 | 34.35 | 35.05 | -118.69 | -128.44 | -209.89 | -29.55 | -18.56 | 0.03 | 0.01 | 0.07 | -0.02 | -0.03 | 0.00 | -0.04 | -0.01 | -0.25 | 0.01 | -0.77 | 0.40 |
| attSPL | -205.08 | -0.66 | -177.56 | -26.86 | -0.08 | 35.05 | 37.49 | 33.06 | 34.36 | -119.69 | 94.87 | -274.74 | 122.55 | -30.48 | 0.04 | -0.03 | 0.10 | -0.06 | -0.05 | 0.00 | -0.06 | -0.01 | 2.08 | 0.38 | -0.61 | 1.09 |
| decSPL | -41.76 | 0.25 | -36.46 | -5.55 | 0.10 | 34.35 | 33.06 | 34.70 | 34.18 | -130.19 | -122.66 | -211.16 | -48.42 | -5.80 | 0.03 | 0.02 | 0.04 | 0.01 | -0.02 | 0.00 | -0.03 | 0.00 | -1.13 | -0.10 | -0.53 | 0.06 |
| earlyDecaySPL | -90.45 | 0.43 | -76.60 | -14.28 | 0.14 | 35.05 | 34.36 | 34.18 | 34.88 | -119.50 | -177.76 | -202.01 | -60.17 | -15.25 | 0.02 | 0.02 | 0.06 | -0.01 | -0.02 | -0.01 | -0.03 | -0.01 | -0.64 | -0.06 | -0.63 | 0.24 |
| totSC | -449.09 | 8.36 | -325.99 | -131.46 | 1.11 | -118.69 | -119.69 | -130.19 | -119.50 | 5907.01 | 1707.28 | 8098.86 | 1471.12 | -94.91 | -0.64 | 0.24 | -0.75 | -0.58 | 0.58 | 0.17 | 0.72 | 0.12 | 9.18 | -3.27 | 1.29 | 5.57 |
| attSC | -7391.47 | -171.34 | -7202.60 | -17.53 | -37.72 | -128.44 | 94.87 | -122.66 | -177.76 | 1707.28 | 62304.01 | -5496.05 | 21851.87 | -301.21 | 5.33 | -5.54 | 6.07 | -4.17 | -1.74 | 1.92 | -2.23 | -0.10 | 314.17 | 49.56 | -7.49 | 83.34 |
| decSC | 1448.04 | 47.75 | 1404.24 | -3.96 | 8.53 | -209.89 | -274.74 | -211.16 | -202.01 | 8098.86 | -5496.05 | 13980.47 | -3473.23 | 115.29 | -1.76 | 1.22 | -2.51 | 0.47 | 1.42 | 0.13 | 1.80 | 0.24 | -53.07 | -16.77 | -0.43 | -13.80 |
| earlyDecaySC | -6574.23 | -101.81 | -5920.35 | -552.06 | -19.80 | -29.55 | 122.55 | -48.42 | -60.17 | 1471.12 | 21851.87 | -3473.23 | 16789.09 | -752.03 | 1.85 | -2.86 | 2.51 | -3.28 | -1.54 | 0.35 | -1.93 | -0.14 | 171.49 | 29.11 | 4.67 | 61.51 |
| TC | 999.35 | 2.68 | 867.57 | 129.09 | 0.40 | -18.56 | -30.48 | -5.80 | -15.25 | -94.91 | -301.21 | 115.29 | -752.03 | 183.60 | 0.06 | 0.02 | -0.16 | 0.32 | 0.14 | 0.00 | 0.17 | 0.02 | -8.94 | -1.05 | 1.79 | -4.55 |
| totFlat | -0.79 | -0.02 | -0.87 | 0.09 | 0.00 | 0.03 | 0.04 | 0.03 | 0.02 | -0.64 | 5.33 | -1.76 | 1.85 | 0.06 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.04 | 0.01 | -0.01 | 0.01 |
| attFlat | 1.26 | 0.02 | 1.19 | 0.05 | 0.00 | 0.01 | -0.03 | 0.02 | 0.02 | 0.24 | -5.54 | 1.22 | -2.86 | 0.02 | 0.00 | 0.01 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -0.05 | -0.01 | 0.00 | -0.01 |
| decFlat | -2.41 | -0.02 | -2.28 | -0.10 | -0.01 | 0.07 | 0.10 | 0.04 | 0.06 | -0.75 | 6.07 | -2.51 | 2.51 | -0.16 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.06 | 0.01 | -0.01 | 0.02 |
| earlyDecayFlat | 2.57 | 0.01 | 2.27 | 0.28 | 0.00 | -0.02 | -0.06 | 0.01 | -0.01 | -0.58 | -4.17 | 0.47 | -3.28 | 0.32 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -0.05 | -0.01 | 0.00 | -0.02 |
| totSpecFlat | 1.05 | 0.01 | 0.94 | 0.11 | 0.00 | -0.03 | -0.05 | -0.02 | -0.02 | 0.58 | -1.74 | 1.42 | -1.54 | 0.14 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -0.02 | 0.00 | 0.00 | -0.01 |
| attSpecFlat | -0.17 | 0.00 | -0.17 | 0.01 | 0.00 | 0.00 | 0.00 | 0.00 | -0.01 | 0.17 | 1.92 | 0.13 | 0.35 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.01 | 0.00 | 0.00 | 0.00 |
| decSpecFlat | 1.31 | 0.01 | 1.16 | 0.13 | 0.00 | -0.04 | -0.06 | -0.03 | -0.03 | 0.72 | -2.23 | 1.80 | -1.93 | 0.17 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -0.02 | 0.00 | 0.00 | -0.01 |
| earlyDecaySpecFlat | 0.12 | 0.00 | 0.11 | 0.02 | 0.00 | -0.01 | -0.01 | 0.00 | -0.01 | 0.12 | -0.10 | 0.24 | -0.14 | 0.02 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| totCrest | -90.78 | -1.26 | -82.96 | -6.57 | -0.26 | -0.25 | 2.08 | -1.13 | -0.64 | 9.18 | 314.17 | -53.07 | 171.49 | -8.94 | 0.04 | -0.05 | 0.06 | -0.05 | -0.02 | 0.01 | -0.02 | 0.00 | 2.99 | 0.50 | -0.04 | 0.88 |
| attCrest | -13.10 | -0.26 | -12.27 | -0.57 | -0.05 | 0.01 | 0.38 | -0.10 | -0.06 | -3.27 | 49.56 | -16.77 | 29.11 | -1.05 | 0.01 | -0.01 | 0.01 | -0.01 | 0.00 | 0.00 | 0.00 | 0.00 | 0.50 | 0.11 | -0.01 | 0.14 |
| decCrest | 27.04 | -0.03 | 25.25 | 1.83 | -0.01 | -0.77 | -0.61 | -0.53 | -0.63 | 1.29 | -7.49 | -0.43 | 4.67 | 1.79 | -0.01 | 0.00 | -0.01 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -0.04 | -0.01 | 0.29 | -0.04 |
| earlyDecayCrest | -39.39 | -0.40 | -35.30 | -3.68 | -0.08 | 0.40 | 1.09 | 0.06 | 0.24 | 5.57 | 83.34 | -13.80 | 61.51 | -4.55 | 0.01 | -0.01 | 0.02 | -0.02 | -0.01 | 0.00 | -0.01 | 0.00 | 0.88 | 0.14 | -0.04 | 0.34 |
# heatmap(cov_mat)
# image(cov_mat)
cov_mat.z <- scale(cov_mat, center = T, scale = T)
cov_mat.z_melt <- reshape2::melt(cov_mat.z)
ggplot(cov_mat.z_melt, aes(Var1, Var2)) +
geom_tile(aes(fill = value)) +
scale_fill_gradient2(low = "red", mid = "orange", high = "yellow") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
axis.title = element_blank())
2 Modeling
Sofia, on 2020-07-07: “The main hypothesis is that playing instruction (N/C) will affect the stroke in a way that is perceivable”.
2.1 Univariate Models –
transSC
Looking at the variables agreed upon I have
decided to use as response variables in the regression
models—implying that this is by no means set in stone—,
it seems like a skewed normal link function would be appropriate to
model them.
Let’s start with
transSCtotDur
(m0_form <-bf(transSC ~ 1 + (1 | subj)))transSC ~ 1 + (1 | subj) which claims that transSC is explained by (‘~’) an
intercept, denoted by ‘1’, and an additional term ‘(1 | subj)’. The ‘1’
in parentheses again stands for the intercept, but the pipe ‘|’ assigns
an intercept to each level of the factor ‘subj’. In this particular case
this means that the model will estimate an individual intercept for each
unique drummer listed in the data set column subj. These
individual, or varying, intercepts are then used in informing the
estimation of the population intercept.
Note: set MODEL to TRUE at the top of the
script if you haven’t compiled/built your model yet.
if (MODEL) {
m0 <- brm(m0_form,
family = skew_normal(),
init = "0",
data = drum_beats)
m0 <- add_criterion(m0,
"loo",
reloo = TRUE)
save(m0,
file = "m0.rda")
} else {
load("m0.rda")
}This null model is also termed unconditional model because it has no grouping structure apart from individuals–the ‘subj’ bit in the model equation above. There is some variation in every natural data set. To make sure, it’s not just variation caused by different participants, we can calculate the intra-class correlation coefficient (ICC).
m0_icc <- ICC(m0, "subj")The Null model’s ICC amounts to 0.76, which suggests that approx. 76 percent of the variation in the data set can be attributed to (or explained by) the grouping structure. This is highly unfortunate, as it does not leave a lot of variation to be explained by independent factors like instruction, or arm. The most likely reason for this high ICC value is the small sample size combined with high inter-individual variation. Small sample sizes combined with large trial numbers are less of a problem when subjects respond, on average, close to the population mean, even with large spread due to fluctuating alertness, increasing fatigue etc., .e.g. in reaction time paradigms. But here, with large intra- and inter-individual variation, this might become a problem.
Tab. 1. Model summary.
(m0_summary <- summary(m0,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ 1 + (1 | subj)
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
alpha ~ normal(0, 4)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ student_t(3, 0, 136.5)
<lower=0> sigma ~ student_t(3, 0, 136.5)
Multilevel Hyperparameters:
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 116.78 34.58 69.10 201.98 1.01 690 699
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 747.16 39.66 665.75 828.48 1.01 703 962
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 64.98 1.50 62.02 68.02 1.00 2270 2035
alpha 3.14 0.33 2.53 3.85 1.00 1849 1521
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The intercept in Tab. 1 amounts to roughly
747. In this model specification, the intercept is identical to the data
set average. The table also shows that the inter-individual standard
deviation (sd(Intercept)) is large compared to the
unexplained variation (sigma). This led to the large ICC
value above. By adding more and more independent factors to the model
specification, we will later try to decrease \(\sigma\), i.e. ‘explain away’ as
much remaining variation as possible.
pp_check(m0, ndraws = 100)
Fig. 12. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
The posterior predictive plot in Fig. 12 gives an impression on how the null model would generate data, given the parameters it estimated from the empirical data. As the thick line deviates from the modeled thin lines, especially on the left and the right side of the peak of the distribution, it is apparent that the model can be improved.
# conditional_effects(m0) 2.1.1 Group Model
In this model, the intercept is complemented with second ‘main’ or
population effect, the grouping variable patient:
(m1_form <- bf(transSC ~ 1 + patient +
(1 | subj)))transSC ~ 1 + patient + (1 | subj) I could have also specified the model without the explicit ‘1 +’, as the intercept is implicitly included in the model unless I explicitly exclude it. From now on I will always save some extra typing by refraining from explicitly indluding the intercept in models.
So now the model not only contains the individuals as grouping structure to ‘explain away’ variation, but also whether they belong to the patients or the healthy subjects.
The prior distributions in Bayesian models reflect the knowledge
about the estimated parameters. The package brms
automatically places weakly informative priors on parameters as soft
contraints. But with more complex models involving varying parameter
estimates, the statistical back end which does the heavy lifting, needs
stronger priors particularly on these parameters, otherwise models don’t
converge. The parameters most vulnerable to outliers in the data are the
varying effects parameters, or random effects in traditional statistics.
Therefore we place a stronger prior probability distribution over the
estimate of the SD of individual intercepts:
(m1_prior <- set_prior("normal(0, 10)", class = "sd"))sd ~ normal(0, 10)and leave the rest of the priors as suggested by the package
brms (see Tab. 2 for their
priors).
if (MODEL) {
m1 <- brm(m1_form,
prior = m1_prior,
init = "0",
family = skew_normal(),
data = drum_beats)
m1 <- add_criterion(m1,
"loo",
reloo = TRUE)
save(m1,
file = "m1.rda")
} else {
load("m1.rda")
}Tab. 2. Model summary.
(m1_summary <- summary(m1,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ 1 + patient + (1 | subj)
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
alpha ~ normal(0, 4)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 10)
<lower=0> sigma ~ student_t(3, 0, 136.5)
Multilevel Hyperparameters:
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 49.11 4.74 40.50 59.08 1.00 2105 2383
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 762.85 24.82 715.40 812.67 1.00 1340 1601
patient1 -26.17 34.41 -93.26 40.92 1.00 1352 1536
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 65.12 1.53 62.26 68.27 1.00 3256 2760
alpha 3.26 0.33 2.64 3.94 1.00 3283 2630
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The intercept in Tab. 2 amounts to roughly
763. In this model specification, the intercept is identical to the
controlled strokes, while the patient1 value is the mean of
the posterior distribution difference between the healthy
subjects and the patients (-26). The table also shows that the
interindividual standard deviation (sd(Intercept) \(\approx\) 49 is still large compared to the
unexplained variation (sigma \(\approx\) 65).
pp_check(m1, ndraws = 100)
Fig. 13. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
The posterior predictive plot in Fig. 13 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 12, which is not very surprising given the small estimated difference between the groups.
conditional_effects(m1,
dpar = "mu")
2.1.2 Group + Condition Model
In this model, the intercept is complemented with second ‘main’ or
population effect, the instruction condition cond:
(m2_form <-bf(transSC ~ patient + cond +
(1 | subj)))transSC ~ patient + cond + (1 | subj) So now the model also contains the individuals as grouping structure to ‘explain away’ variation, but also the manipulations.
The prior distributions in Bayesian models reflect the knowledge
about the estimated parameters. The package brms
automatically places weakly informative priors on parameters as soft
contraints. But with more complex models involving varying parameter
estimates, the statistical back end which does the heavy lifting, needs
stronger priors particularly on these parameters, otherwise models don’t
converge. The parameters most vulnerable to outliers in the data are the
varying effects parameters, or random effects in traditional statistics.
Therefore we place a stronger prior probability distribution over the
estimate of the SD of individual intercepts:
(m2_prior <- c(set_prior("normal(0, 3)", class = "sd"),
set_prior("normal(0, 3)", class = "sigma"),
set_prior("normal(0, 2)", class = "alpha"))) prior class coef group resp dpar nlpar lb ub source
normal(0, 3) sd <NA> <NA> user
normal(0, 3) sigma <NA> <NA> user
normal(0, 2) alpha <NA> <NA> userand leave the rest of the priors as suggested by the package
brms (see Tab. 3 for their
priors).
if (MODEL) {
m2 <- brm(m2_form,
prior = m2_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m2 <- add_criterion(m2,
"loo",
reloo = TRUE)
save(m2,
file = "m2.rda")
} else {
load("m2.rda")
}Tab. 3. Model summary.
(m2_summary <- summary(m2,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient + cond + (1 | subj)
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)
Multilevel Hyperparameters:
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.25 1.51 26.45 32.24 1.00 2829 2544
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 779.96 15.20 750.21 808.91 1.01 1477 2027
patient1 -27.00 20.87 -67.13 14.96 1.01 1390 2000
condN -41.31 3.43 -47.98 -34.63 1.00 4054 3059
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.37 0.98 50.48 54.27 1.00 4593 2495
alpha 1.80 0.24 1.36 2.28 1.00 4153 2874
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The intercept in Tab. 3 amounts to roughly
780. In this model specification, the intercept is identical to the
controlled strokes, while the condN value is the mean of
the posterior distribution difference between controlled and
normal strokes (-41). The table also shows that the interindividual
standard deviation (sd(Intercept) \(\approx\) 29 is still large compared to the
unexplained variation (sigma \(\approx\) 52).
pp_check(m2, ndraws = 100)
Fig. 14. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
The posterior predictive plot in Fig. 14 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 12, which is not very surprising given the small estimated difference between the conditions.
beautify_my_plot(plot(conditional_effects(m2,
dpar = "mu"),
plot = FALSE))
2.1.3 Group + ConditionVI Model
The last model included the instruction condition, more realistically
reflecting the true structure of the data set. But it did not
acknowledge that each subject executed several strokes in each of these
conditions. This model includes a term with a varying intercept (VI) for
condition to reflect just that, which will also assist in more
realistically estimate the population effect of cond:
(m3_form <- bf(transSC ~ patient + cond +
(1 | subj) +
(1 | cond)))transSC ~ patient + cond + (1 | subj) + (1 | cond) Including more varying parameters in the model requires the prior on them to be even stronger:
(m3_prior <- m2_prior) prior class coef group resp dpar nlpar lb ub source
normal(0, 3) sd <NA> <NA> user
normal(0, 3) sigma <NA> <NA> user
normal(0, 2) alpha <NA> <NA> userif (MODEL) {
m3 <- brm(m3_form,
prior = m3_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m3 <- add_criterion(m3,
"loo",
reloo = TRUE)
save(m3,
file = "m3.rda")
} else {
load("m3.rda")
}Tab. 4. Model summary.
(m3_summary <- summary(m3,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient + cond + (1 | subj) + (1 | cond)
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)
Multilevel Hyperparameters:
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.42 1.79 0.10 6.62 1.00 2230 1513
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.27 1.50 26.47 32.18 1.00 2628 2689
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 779.20 14.91 750.33 808.74 1.00 1685 2108
patient1 -25.86 20.09 -65.66 13.54 1.00 1704 2393
condN -41.24 5.40 -52.30 -30.02 1.00 2838 2301
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.36 0.93 50.65 54.26 1.00 4618 2948
alpha 1.80 0.23 1.35 2.26 1.00 5002 2914
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The intercept in Tab. 4 amounts to roughly
779. In this model specification, the intercept is identical to the
controlled strokes, while the condN value is the mean of
the posterior distribution difference between controlled and
normal strokes (-41). The table also shows that the inter-individual
standard deviation (sd(Intercept) \(\approx\) 29 is not large anymore compared
to the unexplained variation (sigma \(\approx\) 52).
pp_check(m3, ndraws = 100)
Fig. 15. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
The posterior predictive plot in Fig. 15 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 12, which is not very surprising given the small estimated difference between the conditions.
beautify_my_plot(plot(conditional_effects(m3,
dpar = "mu"),
plot = FALSE))
2.1.4 Group x ConditionVI
The last model included the instruction condition, more realistically
reflecting the true structure of the data set. But it did not
acknowledge that each subject executed several strokes in each of these
conditions. This model includes a term with a varying intercept (VI) for
condition to reflect just that, which will also assist in more
realistically estimate the population effect of cond:
(m4_form <- bf(transSC ~ patient * cond +
(1 | subj) +
(1 | cond)))transSC ~ patient * cond + (1 | subj) + (1 | cond) Including more varying parameters in the model requires the prior on them to be even stronger:
(m4_prior <- m3_prior) prior class coef group resp dpar nlpar lb ub source
normal(0, 3) sd <NA> <NA> user
normal(0, 3) sigma <NA> <NA> user
normal(0, 2) alpha <NA> <NA> userif (MODEL) {
m4 <- brm(m4_form,
prior = m4_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m4 <- add_criterion(m4,
"loo",
reloo = TRUE)
save(m4,
file = "m4.rda")
} else {
load("m4.rda")
}Tab. 5. Model summary.
(m4_summary <- summary(m4,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient * cond + (1 | subj) + (1 | cond)
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)
Multilevel Hyperparameters:
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.43 1.84 0.08 6.90 1.00 2489 1546
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.23 1.48 26.36 32.17 1.00 3437 3199
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 780.51 15.43 751.24 811.37 1.00 1548 2068
patient1 -28.11 21.00 -70.61 11.97 1.00 1554 2272
condN -42.62 6.34 -55.38 -29.58 1.00 2787 2530
patient1:condN 2.44 6.33 -9.93 14.50 1.00 5037 3014
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.39 0.91 50.66 54.21 1.00 5152 3046
alpha 1.81 0.24 1.36 2.31 1.00 4824 3033
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The intercept in Tab. 5 amounts to roughly
781. In this model specification, the intercept is identical to the
controlled strokes, while the condN value is the mean of
the posterior distribution difference between controlled and
normal strokes (-43). The table also shows that the inter-individual
standard deviation (sd(Intercept) \(\approx\) 29 is not large anymore compared
to the unexplained variation (sigma \(\approx\) 52).
pp_check(m4, ndraws = 100)
Fig. 16. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
The posterior predictive plot in Fig. 16 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 12, which is not very surprising given the small estimated difference between the conditions.
beautify_my_plot(plot(conditional_effects(m4,
dpar = "mu"),
plot = FALSE))
2.1.5 Group x ConditionVI + Arm Model
(m5_form <- bf(transSC ~ patient * cond + arm +
(1 | subj) +
(1 | cond)))transSC ~ patient * cond + arm + (1 | subj) + (1 | cond) (m5_prior <- m4_prior) prior class coef group resp dpar nlpar lb ub source
normal(0, 3) sd <NA> <NA> user
normal(0, 3) sigma <NA> <NA> user
normal(0, 2) alpha <NA> <NA> userif (MODEL) {
m5 <- brm(m5_form,
prior = m5_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m5 <- add_criterion(m5,
"loo",
reloo = TRUE)
save(m5,
file = "m5.rda")
} else {
load("m5.rda")
}Tab. 6. Model summary.
(m5_summary <- summary(m5,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient * cond + arm + (1 | subj) + (1 | cond)
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)
Multilevel Hyperparameters:
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.38 1.81 0.10 6.85 1.00 2808 1452
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.26 1.48 26.54 32.26 1.00 3002 2712
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 778.69 15.22 750.01 808.84 1.00 1599 2253
patient1 -27.48 21.16 -67.63 13.57 1.00 1482 1997
condN -42.50 6.19 -54.84 -30.18 1.00 2940 2333
armND 1.94 3.17 -4.19 8.16 1.00 4885 2750
patient1:condN 2.44 6.31 -9.70 15.13 1.00 3896 2989
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.39 0.93 50.62 54.24 1.00 4932 3065
alpha 1.80 0.23 1.36 2.26 1.00 4264 2854
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Population-Level Effects section of Tab.
6 are now two estimates for the skewness parameter alpha (one for
the controlled [alpha_Intercept] and one for the normal
condition [alpha_condN]), while there is no entry for alpha
in the Family Specific section anymore. All this due to our explicit
modeling of alpha conditional on instruction condition.
pp_check(m5, ndraws = 100)
Fig. 17. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
beautify_my_plot(plot(conditional_effects(m5,
dpar = "mu"),
plot = FALSE))
2.1.6 Group x Condition + ArmVI Model
(m6_form <-bf(transSC ~ patient * cond + arm +
(1 | subj) +
(1 | cond) +
(1 | arm)))transSC ~ patient * cond + arm + (1 | subj) + (1 | cond) + (1 | arm) (m6_prior <- m5_prior) prior class coef group resp dpar nlpar lb ub source
normal(0, 3) sd <NA> <NA> user
normal(0, 3) sigma <NA> <NA> user
normal(0, 2) alpha <NA> <NA> userif (MODEL) {
m6 <- brm(m6_form,
prior = m6_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m6 <- add_criterion(m6,
"loo",
reloo = TRUE)
save(m6,
file = "m6.rda")
} else {
load("m6.rda")
}Tab. 7. Model summary.
(m6_summary <- summary(m6,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient * cond + arm + (1 | subj) + (1 | cond) + (1 | arm)
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)
Multilevel Hyperparameters:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.36 1.80 0.09 6.72 1.00 1732 1395
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.37 1.78 0.09 6.61 1.00 2597 1930
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.27 1.51 26.41 32.38 1.00 3277 2942
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 778.89 16.31 746.65 810.19 1.00 1636 1936
patient1 -27.71 21.59 -69.78 14.20 1.00 1737 1690
condN -42.73 6.05 -55.41 -30.77 1.00 2625 2270
armND 1.96 5.20 -8.77 12.69 1.00 2936 2231
patient1:condN 2.57 6.17 -9.32 14.65 1.00 4105 2772
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.36 0.93 50.50 54.21 1.00 4749 2764
alpha 1.80 0.23 1.35 2.27 1.00 5175 3159
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Population-Level Effects section of Tab.
7 are now two estimates for the skewness parameter alpha (one for
the controlled [alpha_Intercept] and one for the normal
condition [alpha_condN]), while there is no entry for alpha
in the Family Specific section anymore. All this due to our explicit
modeling of alpha conditional on instruction condition.
pp_check(m6, ndraws = 100)
Fig. 18. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
beautify_my_plot(plot(conditional_effects(m6,
dpar = "mu"),
plot = FALSE))
2.1.7 Group x Condition x ArmVI Model
(m6a_form <-bf(transSC ~ patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm)))transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) (m6a_prior <- m6_prior) prior class coef group resp dpar nlpar lb ub source
normal(0, 3) sd <NA> <NA> user
normal(0, 3) sigma <NA> <NA> user
normal(0, 2) alpha <NA> <NA> userif (MODEL) {
m6a <- brm(m6a_form,
prior = m6a_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m6a <- add_criterion(m6a,
"loo",
reloo = TRUE)
save(m6a,
file = "m6a.rda")
} else {
load("m6a.rda")
}Tab. 8. Model summary.
(m6a_summary <- summary(m6a,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)
Multilevel Hyperparameters:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.38 1.82 0.07 6.83 1.00 2462 1425
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.44 1.84 0.10 6.72 1.00 2523 1781
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 29.39 1.47 26.64 32.40 1.00 3165 2828
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 784.68 15.62 754.83 815.57 1.00 1612
patient1 -42.15 21.48 -83.41 -0.29 1.00 1582
condN -42.56 7.49 -57.33 -27.57 1.00 2083
armND -7.65 7.22 -22.05 6.73 1.00 2220
patient1:condN 9.10 8.96 -8.29 26.78 1.00 2454
patient1:armND 26.71 8.44 10.13 43.50 1.00 2613
condN:armND -0.19 8.56 -17.05 16.61 1.00 2081
patient1:condN:armND -13.95 11.97 -38.00 9.13 1.00 2073
Tail_ESS
Intercept 2287
patient1 2147
condN 2445
armND 2817
patient1:condN 2823
patient1:armND 2724
condN:armND 2508
patient1:condN:armND 2766
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 52.16 0.91 50.44 54.03 1.00 4903 3201
alpha 1.68 0.24 1.21 2.18 1.00 4802 2979
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Population-Level Effects section of
Tab. 8 are now two estimates for the skewness
parameter alpha (one for the controlled [alpha_Intercept]
and one for the normal condition [alpha_condN]), while
there is no entry for alpha in the Family Specific section anymore. All
this due to our explicit modeling of alpha conditional on instruction
condition.
pp_check(m6a, ndraws = 100)
Fig. 19. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
beautify_my_plot(plot(conditional_effects(m6a,
dpar = "mu"),
plot = FALSE))
conditions <- make_conditions(drum_beats,
vars = "arm")
beautify_my_plot(plot(conditional_effects(m6a,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
2.1.8 Modeling Condition-specific Variation
As is apparent from Fig. 1, the controlled condition yielded broader density distributions in most descriptors. Hence, we model condition-dependent variation in the next model:
(m7_form <-bf(transSC ~ patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm),
sigma ~ cond))transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ condThis model has now two outcomes, not just one, as the models before.
Spread is estimated as sigma, and it varies conditional on
instruction condition.
Here’s the non-standard (additional) prior:
(m7_prior <- c(set_prior("normal(0, 5)", class = "sd"),
set_prior("normal(0, 2)", class = "alpha"),
set_prior("normal(0, 10)", class = "Intercept", dpar = "sigma"),
set_prior("normal(0, 10)", class = "b", dpar = "sigma")
)) prior class coef group resp dpar nlpar lb ub source
normal(0, 5) sd <NA> <NA> user
normal(0, 2) alpha <NA> <NA> user
normal(0, 10) Intercept sigma <NA> <NA> user
normal(0, 10) b sigma <NA> <NA> userif (MODEL) {
m7 <- brm(m7_form,
prior = m7_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m7 <- add_criterion(m7,
"loo",
reloo = TRUE)
save(m7,
file = "m7.rda")
} else {
load("m7.rda")
}Tab. 9. Model summary.
(m7_summary <- summary(m7,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
alpha ~ normal(0, 2)
b_sigma ~ normal(0, 10)
Intercept ~ student_t(3, 735, 136.5)
Intercept_sigma ~ normal(0, 10)
<lower=0> sd ~ normal(0, 5)
Multilevel Hyperparameters:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 4.01 3.05 0.13 11.09 1.00 2457 1648
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 3.97 2.97 0.15 11.22 1.00 3501 2247
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 35.46 2.38 30.98 40.35 1.00 4063 3041
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 788.91 20.56 747.69 830.07 1.00 1641
sigma_Intercept 4.30 0.04 4.23 4.38 1.00 5439
patient1 -41.29 26.62 -93.30 11.66 1.00 1621
condN -45.64 9.78 -64.94 -26.17 1.00 2365
armND -9.12 10.98 -30.01 13.56 1.00 2578
patient1:condN 7.46 10.14 -11.97 27.61 1.00 3092
patient1:armND 24.81 11.66 2.43 48.21 1.00 2642
condN:armND 0.11 9.70 -18.49 18.96 1.00 2610
patient1:condN:armND -10.28 13.92 -37.46 17.60 1.00 2386
sigma_condN -0.46 0.06 -0.57 -0.34 1.00 4995
Tail_ESS
Intercept 2320
sigma_Intercept 2463
patient1 2120
condN 2489
armND 1771
patient1:condN 2831
patient1:armND 2941
condN:armND 2884
patient1:condN:armND 2882
sigma_condN 2696
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
alpha 2.25 0.29 1.72 2.86 1.00 4782 3109
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Links section of the summary (Tab.
9) we note that sigma is no longer modeled on the identity scale but
on the \(\log_{2}\) scale. In the
Population-Level Effects section of Tab. 9 we
now find not only the estimates for the Intercept and condN
but also two estimates for sigma (one for the controlled
[sigma_Intercept] and one for the normal condition
[sigma_condN]). Because of this, there is no entry for
sigma in the Family Specific section anymore. All this due to our
explicit modeling of sigma conditional on instruction condition.
pp_check(m7, ndraws = 100)
Fig. 20. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
beautify_my_plot(plot(conditional_effects(m7,
dpar = "mu"),
plot = FALSE))
conditions <- make_conditions(drum_beats,
vars = "arm")
beautify_my_plot(plot(conditional_effects(m7,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m7,
dpar = "sigma"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m7,
dpar = "sigma",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
Fig. 21. Conditional plot.
In Fig. 21 we see a difference between the
mean (‘mu’) estimated value for transSC depending on
instruction condition (A), but no difference between dominant and
non-dominant arm (B). Consequently, there is no interaction in (C).
The spread (sigma) of the estimated distributions also
differs between conditions (D), but again not between arms (E). The
latter is not surprising because we modeled sigma to vary conditional on
condition, not arm.
2.1.9 Modeling Arm-specific Skewness
At least in some descriptors in Fig. 1 the skewness of the distribution seems to change depending on the arm. The skew normal distribution is a generalization of the Gaussian distribution, allowing for the additional shape-parameter skewness (asymmetry) to vary. Hence, we model this side-dependent skewness in the next model:
(m8_form <-bf(transSC ~ patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm),
sigma ~ cond,
alpha ~ arm))transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond
alpha ~ arm(m8_prior <- c(set_prior("normal(0, 3)", class = "sd"),
set_prior("normal(0, 10)", class = "Intercept", dpar = "sigma"),
set_prior("normal(0, 10)", class = "b", dpar = "sigma"),
set_prior("normal(0, 2)", class = "Intercept", dpar = "alpha"),
set_prior("normal(0, 2)", class = "b", dpar = "alpha"))
) prior class coef group resp dpar nlpar lb ub source
normal(0, 3) sd <NA> <NA> user
normal(0, 10) Intercept sigma <NA> <NA> user
normal(0, 10) b sigma <NA> <NA> user
normal(0, 2) Intercept alpha <NA> <NA> user
normal(0, 2) b alpha <NA> <NA> userif (MODEL) {
m8 <- brm(m8_form,
prior = m8_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m8 <- add_criterion(m8,
"loo",
reloo = TRUE)
save(m8,
file = "m8.rda")
} else {
load("m8.rda")
}Tab. 10. Model summary.
(m8_summary <- summary(m8,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond
alpha ~ arm
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 10)
Intercept ~ student_t(3, 735, 136.5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 10)
<lower=0> sd ~ normal(0, 3)
Multilevel Hyperparameters:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.35 1.78 0.09 6.59 1.00 2975 1680
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.42 1.83 0.10 6.75 1.00 3645 2387
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 27.38 1.51 24.49 30.49 1.00 2970 2278
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 787.52 15.14 757.90 816.85 1.00 1551
sigma_Intercept 4.33 0.04 4.26 4.41 1.00 4488
alpha_Intercept 1.69 0.41 0.84 2.40 1.00 4289
patient1 -43.17 20.57 -82.05 -1.86 1.00 1674
condN -44.43 8.36 -61.25 -28.21 1.00 2119
armND -9.63 9.07 -26.73 8.35 1.00 2117
patient1:condN 9.38 10.65 -11.30 30.57 1.00 2068
patient1:armND 27.49 11.77 3.94 50.13 1.00 1892
condN:armND -1.48 9.57 -20.15 17.16 1.00 1962
patient1:condN:armND -9.96 13.91 -37.79 17.02 1.00 1881
sigma_condN -0.51 0.06 -0.62 -0.40 1.00 4226
alpha_armND 1.68 0.54 0.67 2.80 1.00 3751
Tail_ESS
Intercept 2375
sigma_Intercept 2913
alpha_Intercept 2079
patient1 2309
condN 2255
armND 2509
patient1:condN 2545
patient1:armND 2594
condN:armND 2428
patient1:condN:armND 2228
sigma_condN 3109
alpha_armND 2446
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Population-Level Effects section of Tab.
10 are now two estimates for the skewness parameter alpha (one for
the dominant [alpha_Intercept] and one for the non-dominant
arm [alpha_armND]), while there is no entry for alpha in
the Family Specific section anymore. All this due to our explicit
modeling of alpha conditional on arm.
pp_check(m8, ndraws = 100)
Fig. 22. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
beautify_my_plot(plot(conditional_effects(m8,
dpar = "mu"),
plot = FALSE))
conditions <- make_conditions(drum_beats,
vars = "arm")
beautify_my_plot(plot(conditional_effects(m8,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m8,
dpar = "sigma"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m8,
dpar = "sigma",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m8,
dpar = "alpha"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m8,
dpar = "alpha",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
Fig. 23. Conditional plots.
In Fig. 23 there is a clear difference in estimated mean (‘mu’) of \(transSC\) conditional on instruction (A), but not on arm (B). Consequently, there is no interaction between instruction and arm (C). On the other hand, there is a difference in estimated skewness (‘alpha’) between sides (E), but not instructions (D). (F) follows from that. ‘sigma’ shows a clear difference between instructions (G), not so for side (H), which is mirrored in (I).
2.1.10 Modeling Influence of Arm on Shape
To be able to see whether skewness and spread are depending on both experimental manipulation and arm, the model formulas for ‘sigma’ and ‘alpha’ are updated accordingly:
(m9_form <-bf(transSC ~ patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm),
sigma ~ cond,
alpha ~ arm))transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond
alpha ~ arm(m9_prior <- m8_prior) prior class coef group resp dpar nlpar lb ub source
normal(0, 3) sd <NA> <NA> user
normal(0, 10) Intercept sigma <NA> <NA> user
normal(0, 10) b sigma <NA> <NA> user
normal(0, 2) Intercept alpha <NA> <NA> user
normal(0, 2) b alpha <NA> <NA> userif (MODEL) {
m9 <- brm(m9_form,
prior = m9_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m9 <- add_criterion(m9,
"loo",
reloo = TRUE)
save(m9,
file = "m9.rda")
} else {
load("m9.rda")
}Tab. 11. Model summary.
(m9_summary <- summary(m9,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond
alpha ~ arm
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 10)
Intercept ~ student_t(3, 735, 136.5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 10)
<lower=0> sd ~ normal(0, 3)
Multilevel Hyperparameters:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.41 1.80 0.10 6.62 1.00 3102 1803
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.40 1.81 0.12 6.79 1.00 3362 2082
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 27.40 1.51 24.51 30.42 1.00 3438 3097
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 787.88 15.43 757.82 819.11 1.00 1707
sigma_Intercept 4.33 0.04 4.26 4.41 1.00 4444
alpha_Intercept 1.70 0.40 0.89 2.45 1.00 4143
patient1 -44.33 21.11 -86.36 -2.99 1.00 1676
condN -44.87 8.20 -61.12 -28.53 1.00 2114
armND -10.20 9.12 -27.85 8.26 1.00 1934
patient1:condN 10.29 10.46 -10.00 31.03 1.00 1816
patient1:armND 28.35 11.81 5.73 51.76 1.00 1713
condN:armND -0.81 9.48 -19.85 17.70 1.00 1904
patient1:condN:armND -11.21 13.77 -38.63 15.43 1.00 1541
sigma_condN -0.52 0.06 -0.63 -0.40 1.00 4153
alpha_armND 1.66 0.55 0.60 2.81 1.00 3858
Tail_ESS
Intercept 2215
sigma_Intercept 2965
alpha_Intercept 2447
patient1 2169
condN 2946
armND 2731
patient1:condN 2771
patient1:armND 2531
condN:armND 2694
patient1:condN:armND 2299
sigma_condN 3089
alpha_armND 2679
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the Population-Level Effects section of Tab. 11 are now three estimates for the skewness parameter alpha, as well as three estimates for spread (‘sigma’).
pp_check(m9, ndraws = 100)
Fig. 24. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
beautify_my_plot(plot(conditional_effects(m9,
dpar = "mu"),
plot = FALSE))
conditions <- make_conditions(drum_beats,
vars = "arm")
beautify_my_plot(plot(conditional_effects(m9,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m9,
dpar = "sigma"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m9,
dpar = "sigma",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m9,
dpar = "alpha"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m9,
dpar = "alpha",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
Fig. 25. Conditional plots. Had to take out all the fancy figure polishing, subplot lettering etc. Subplots are supposed to be numbered from A to … starting in the top left corner, proceeding row-wise.
In Fig. 25 the dots are the means of the posterior distributions of the respective estimate, and can be interpreted similarly to the point estimates in traditional frequentist statistics. The error bars represent 95% credible intervals, which are interpreted as comprising the true value with 95% probability, given the model and the data.
There is a difference in estimated the mean (‘mu’) of \(transSC\) conditional on group (A) and condition (B), not so on arm (C). Likewise, there is a two-way interaction between group and arm (D), but not between condition and arm (E). Consequently, there is no three-way interaction (F).
The remaining, unexplained variation ‘sigma’ shows a clear difference between instructions (H) and lesser so for side (I), but albeit their opposing trends, there is no interaction (J thru M).
There is a difference in estimated skewness (‘alpha’) both between instruction (O) and side (P), but with opposing trends. Nevertheless, there is no interaction between them (Q thru T).
From a data-driven perspective this should lead to a simplification of the model, leaving out all the unnecessary interaction terms. A theory-driven approach though requires the final model to be ‘full’ in the sense that all explanatory variables must interact with all others, implying that not finding an interaction between factors either is due to there really not being an interaction–rendering the model mis-specified and proving the underlying theory wrong–, or due to insufficient data. Since the data set at hand is very small I vote to blame the failure to find interactions between all factors on the small sample size and would therefore keep all interactions in the model.
2.1.11 Modeling Interactions of Arm and Instruction on Shape and Spread
If an interaction between cond and arm
should be part of the model explaining the outcome transSC,
then this interactin should be also present in the part of the model
explaining sigma and alpha:
(m10_form <-bf(transSC ~ patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm),
sigma ~ cond * arm,
alpha ~ cond * arm))transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond * arm
alpha ~ cond * arm(m10_prior <- m9_prior) prior class coef group resp dpar nlpar lb ub source
normal(0, 3) sd <NA> <NA> user
normal(0, 10) Intercept sigma <NA> <NA> user
normal(0, 10) b sigma <NA> <NA> user
normal(0, 2) Intercept alpha <NA> <NA> user
normal(0, 2) b alpha <NA> <NA> userif (MODEL) {
m10 <- brm(m10_form,
prior = m10_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m10 <- add_criterion(m10,
"loo",
reloo = TRUE)
save(m10,
file = "m10.rda")
} else {
load("m10.rda")
}Tab. 12. Model summary.
(m10_summary <- summary(m10,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ cond * arm
alpha ~ cond * arm
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 10)
Intercept ~ student_t(3, 735, 136.5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 10)
<lower=0> sd ~ normal(0, 3)
Multilevel Hyperparameters:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.41 1.83 0.13 6.86 1.00 3683 2088
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.44 1.79 0.10 6.69 1.00 2636 1699
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 27.40 1.50 24.65 30.46 1.00 3648 2642
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 787.53 15.38 757.35 817.04 1.00 1984
sigma_Intercept 4.29 0.05 4.20 4.40 1.00 3976
alpha_Intercept 4.40 0.65 3.23 5.72 1.00 4080
patient1 -38.10 20.72 -78.43 2.13 1.00 2100
condN -44.27 7.81 -59.75 -28.73 1.00 2882
armND -9.11 9.43 -27.43 9.42 1.00 2618
patient1:condN 2.15 9.19 -16.33 19.86 1.00 2719
patient1:armND 31.16 10.83 9.38 52.14 1.00 2578
condN:armND -2.98 9.90 -22.31 16.33 1.00 2465
patient1:condN:armND -14.15 13.04 -39.40 11.65 1.00 2271
sigma_condN -0.52 0.07 -0.66 -0.38 1.00 3888
sigma_armND 0.13 0.07 -0.00 0.27 1.00 3423
sigma_condN:armND -0.18 0.09 -0.36 -0.00 1.00 3469
alpha_condN -5.35 0.87 -6.92 -3.46 1.00 2477
alpha_armND 2.05 1.03 0.16 4.16 1.00 3322
alpha_condN:armND 1.05 1.16 -1.35 3.24 1.00 2460
Tail_ESS
Intercept 2510
sigma_Intercept 3305
alpha_Intercept 2860
patient1 2414
condN 2815
armND 3087
patient1:condN 2812
patient1:armND 2815
condN:armND 2630
patient1:condN:armND 2618
sigma_condN 3186
sigma_armND 2930
sigma_condN:armND 2914
alpha_condN 2856
alpha_armND 2450
alpha_condN:armND 2580
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(m10, ndraws = 100)
Fig. 26. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
beautify_my_plot(plot(conditional_effects(m10,
dpar = "mu"),
plot = FALSE))
conditions <- make_conditions(drum_beats,
vars = "arm")
beautify_my_plot(plot(conditional_effects(m10,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m10,
dpar = "sigma"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m10,
dpar = "sigma",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m10,
dpar = "alpha"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m10,
dpar = "alpha",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
Fig. 27. Conditional plots. Had to
take out all the fancy figure polishing. Subplots are supposed to be
numbered from A to … starting in the top left corner, proceeding
row-wise. Instead, they are now numbered from A to F, plus one
un-numbered plot with the three-way interaction, for each of the
estimated outcomes transSC, sigma, and
alpha
In Fig. 27 the dots are the means of the posterior distributions of the respective estimate, and can be interpreted similarly to the point estimates in traditional frequentist statistics. The error bars represent 95% credible intervals, which are interpreted as comprising the true value with 95% probability, given the model and the data.
There is a difference in estimated the mean (‘mu’) of \(transSC\) conditional on group (A) and condition (B), not so on arm (C). Likewise, there is a two-way interaction between group and arm (D), but not between condition and arm (E). Consequently, there is no three-way interaction (F).
The remaining, unexplained variation ‘sigma’ shows a clear difference between instructions (H) and lesser so for side (I), but albeit their opposing trends, there is no interaction (J thru M).
There is a difference in estimated skewness (‘alpha’) both between instruction (O) and side (P), but with opposing trends. Nevertheless, there is no interaction between them (Q thru T).
From a data-driven perspective this should lead to a simplification of the model, leaving out all the unnecessary interaction terms. A theory-driven approach though requires the final model to be ‘full’ in the sense that all explanatory variables must interact with all others, implying that not finding an interaction between factors either is due to there really not being an interaction–rendering the model mis-specified and proving the underlying theory wrong–, or due to insufficient data. Since the data set at hand is very small I vote to blame the failure to find interactions between all factors on the small sample size and would therefore keep all interactions in the model.
2.1.12 The Full Model: Modeling the Interactions of Group, Instruction, and Side on Shape and Spread
(m11_form <-bf(transSC ~ patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm),
sigma ~ patient * cond * arm,
alpha ~ patient * cond * arm))transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ patient * cond * arm
alpha ~ patient * cond * arm(m11_prior <- m10_prior) prior class coef group resp dpar nlpar lb ub source
normal(0, 3) sd <NA> <NA> user
normal(0, 10) Intercept sigma <NA> <NA> user
normal(0, 10) b sigma <NA> <NA> user
normal(0, 2) Intercept alpha <NA> <NA> user
normal(0, 2) b alpha <NA> <NA> userif (MODEL) {
m11 <- brm(m11_form,
prior = m11_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m11 <- add_criterion(m11,
"loo",
reloo = TRUE)
save(m11,
file = "m11.rda")
} else {
load("m11.rda")
}Tab. 13. Model summary.
(m11_summary <- summary(m11,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ patient * cond * arm
alpha ~ patient * cond * arm
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 10)
Intercept ~ student_t(3, 735, 136.5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 10)
<lower=0> sd ~ normal(0, 3)
Multilevel Hyperparameters:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.46 1.85 0.12 6.94 1.00 4068 2366
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.40 1.81 0.12 6.90 1.00 4364 2706
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 27.31 1.49 24.51 30.26 1.00 6072 3196
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 790.58 15.68 759.87 820.68 1.00 1928
sigma_Intercept 4.34 0.07 4.21 4.47 1.00 2907
alpha_Intercept 4.42 0.77 3.00 6.00 1.00 4572
patient1 -46.81 21.23 -88.55 -5.94 1.00 2013
condN -47.05 8.30 -63.05 -30.07 1.00 2747
armND -12.06 9.99 -31.54 7.96 1.00 2188
patient1:condN 10.77 9.98 -9.04 30.28 1.00 2459
patient1:armND 41.45 12.83 16.20 66.79 1.00 2070
condN:armND 1.32 10.43 -19.04 21.39 1.00 2085
patient1:condN:armND -26.89 14.85 -55.77 2.18 1.00 1957
sigma_patient1 -0.15 0.10 -0.35 0.05 1.00 2531
sigma_condN -0.54 0.11 -0.74 -0.32 1.00 3088
sigma_armND 0.05 0.10 -0.16 0.23 1.00 2477
sigma_patient1:condN 0.11 0.15 -0.18 0.40 1.00 2678
sigma_patient1:armND 0.21 0.14 -0.06 0.49 1.00 2243
sigma_condN:armND -0.03 0.13 -0.29 0.22 1.00 2595
sigma_patient1:condN:armND -0.33 0.19 -0.69 0.04 1.00 2537
alpha_patient1 -0.61 1.03 -2.64 1.46 1.00 4839
alpha_condN -5.07 1.10 -7.15 -2.98 1.00 2807
alpha_armND 0.64 1.20 -1.60 3.03 1.00 4209
alpha_patient1:condN 0.49 1.21 -1.86 2.88 1.00 3406
alpha_patient1:armND 2.95 1.32 0.41 5.53 1.00 4069
alpha_condN:armND 1.25 1.31 -1.35 3.80 1.00 3592
alpha_patient1:condN:armND -1.36 1.37 -4.08 1.36 1.00 4208
Tail_ESS
Intercept 2519
sigma_Intercept 3364
alpha_Intercept 3149
patient1 2584
condN 2852
armND 2727
patient1:condN 3169
patient1:armND 2793
condN:armND 2875
patient1:condN:armND 2639
sigma_patient1 3017
sigma_condN 3142
sigma_armND 2916
sigma_patient1:condN 2823
sigma_patient1:armND 2913
sigma_condN:armND 3185
sigma_patient1:condN:armND 2737
alpha_patient1 2941
alpha_condN 3011
alpha_armND 3205
alpha_patient1:condN 3008
alpha_patient1:armND 3298
alpha_condN:armND 3017
alpha_patient1:condN:armND 3177
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(m11, ndraws = 100)
Fig. 28. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
beautify_my_plot(plot(conditional_effects(m11,
dpar = "mu"),
plot = FALSE))
conditions <- make_conditions(drum_beats,
vars = "arm")
beautify_my_plot(plot(conditional_effects(m11,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11,
dpar = "sigma"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11,
dpar = "sigma",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11,
dpar = "alpha"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11,
dpar = "alpha",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
Fig. 29. Conditional plots. Had to take out all the fancy figure polishing, subplot lettering etc. Subplots are supposed to be numbered from A to … starting in the top left corner, proceeding row-wise.
In Fig. 29 the dots are the means of the posterior distributions of the respective estimate, and can be interpreted similarly to the point estimates in traditional frequentist statistics. The error bars represent 95% credible intervals, which are interpreted as comprising the true value with 95% probability, given the model and the data.
There is a difference in estimated the mean (‘mu’) of \(transSC\) conditional on group (A) and condition (B), not so on arm (C). Likewise, there is a two-way interaction between group and arm (D), but not between condition and arm (E). Consequently, there is no three-way interaction (F).
The remaining, unexplained variation ‘sigma’ shows a clear difference between instructions (H) and lesser so for side (I), but albeit their opposing trends, there is no interaction (J thru M).
There is a difference in estimated skewness (‘alpha’) both between instruction (O) and side (P), but with opposing trends. Nevertheless, there is no interaction between them (Q thru T).
From a data-driven perspective this should lead to a simplification of the model, leaving out all the unnecessary interaction terms. A theory-driven approach though requires the final model to be ‘full’ in the sense that all explanatory variables must interact with all others, implying that not finding an interaction between factors either is due to there really not being an interaction–rendering the model mis-specified and proving the underlying theory wrong–, or due to insufficient data. Since the data set at hand is very small I vote to blame the failure to find interactions between all factors on the small sample size and would therefore keep all interactions in the model.
2.1.13 Model Comparison
We compare models by their estimated log-posterior density (elpd). The smaller this value the better a model predicts the data, despite the penalty for additional covariates. When the difference between two models is more than two SE apart they are considered to be ‘different enough’ to warrant acceptance of one over the other despite possibly smaller parsimony.
univ_model_compar <- loo_compare(m0, m2, m3, m5, m6, m6a, m7, m8, m9, m10, m11)
print(univ_model_compar, simplify = T) elpd_diff se_diff
m10 0.0 0.0
m11 -2.1 3.7
m9 -40.7 9.0
m8 -40.8 8.9
m7 -43.6 10.0
m6a -83.1 15.8
m3 -83.2 16.7
m2 -83.5 16.7
m6 -85.1 16.8
m5 -85.4 16.8
m0 -121.7 15.5 The model with the smallest elpd is m10,
which corresponds to the second to last model (see Modeling
Interactions of Arm and Instruction on Shape and Spread). The
difference between the best and the second-best model (m11, [Modeling
the Influence of Group on Shape and Spread]) is -2.1, which is less than
roughly two SE (2 x 3.67 = 7.34) away from 0. Therefore, we consider the
best fitting and the second-best fitting models (roughly) equivalent.
Based on theoretical grounds, modeling the three estimates
mu, sigma and alpha should be
done with the full model–i.e., containing all the main effect
interactions–if there is no risk of overfitting. Therefore, in the
following we use the full model ([Modeling the Influence of Group on
Shape and Spread]).
2.1.14 Conclusion
The fact that the full model turned out to be the second-best suggests (see above) that the data really are best fit by this model, and that the patient term and the interaction terms really should remain in the model.
wrap_plots(plot(conditional_effects(m11,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
# facet_args = list(arm = "label_value"),
plot = FALSE)[[1]] +
scale_color_colorblind() +
theme(legend.position = "none"),
plot(conditional_effects(m11,
dpar = "sigma",
effects = "patient:cond",
conditions = conditions),
plot = FALSE)[[1]] +
scale_color_colorblind() +
theme(legend.position = "none"),
plot(conditional_effects(m11,
dpar = "alpha",
effects = "patient:cond",
conditions = conditions),
plot = FALSE)[[1]]) +
scale_color_colorblind() +
plot_annotation(tag_levels = "A")
Fig. 30. Three-way interaction plots of the three modeled parameters \(\mu\), \(\sigma\), and \(\alpha\) as estimated by the full model. The x axis shows group membership (0 = healthy subjects, 1 = patients). The instruction conditions are color-coded (black = controlled strokes, brown = normal strokes), and the facets show the side (1 = dominant, 2 = non-dominant).
Fig. 30A suggests that \(\mu\) is higher for healthy subjects than
for patients, and higher for controlled strokes than for normal strokes
which is really odd on an intuitive level and seems contradictive (at
least to me). I would’ve expected transSC to be at some
optimal level in healthy subjects while performing normal strokes with
their dominant side, and decrease from there. Yet it does not seem to be
so simple. The biggest and clearest difference between any of the
experimental factors is between instruction conditions “normal” and
“controlled” (see Tab. 13 for the credible
intervals of these estimates). Of the interactions shown in
Fig. 30A, only the estimate for
patient1:armND is sufficiently far from zero, while the
central 95% of the posterior distribution of
patient1:condN:armND is, after all, located mostly to the
left of zero.
Fig. 30B might clarify why this is the case, because it shows a higher unexplained variation in the data for controlled strokes as compared to normal strokes which was expected due to controlled strokes are executed with less automatized motor programs (at least that’s what the lay person at the keyboard thinks). And since higher variation is often accompanied by more extreme values, these extreme values might drive the mean of the posterior distribution of \(\mu\) up. Interestingly, \(\sigma\) is estimated to have almost the same low value across groups and sides in the normal instruction condition, whereas it is unanimously high in the controlled condition.
Fig. 30C gives the
impression that in the normal condition the distribution of
transSC is roughly symmetrical around zero, while in the
controlled condition there is an increased amount of higher
transSC values making the distribution heavy-tailed on the
right side (see e.g., Wikipedia on
skewness).
Summarizing these results, transSC is, on average, lower
in normal strokes while approximately normally distributed, and with
lower spread than in controlled strokes. The fact that there is a
difference in the location parameter \(\mu\) between the two instructions suggests
that this descriptor might serve as a cue in listening tests.
2.2 Univariate Models –
transFlat
2.2.1 The Full Model: Modeling the Interactions of Group, Instruction, and Side on Shape and Spread
(m11_transFlat_form <-bf(transFlat ~ 0 + Intercept + patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm),
sigma ~ patient * cond * arm,
alpha ~ patient * cond * arm))transFlat ~ 0 + Intercept + patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ patient * cond * arm
alpha ~ patient * cond * arm(m11_transFlat_prior <- c(set_prior("normal(0.8, 10)", class = "b", coef = "Intercept"),
set_prior("normal(0, 10)", class = "b"),
set_prior("normal(0, 5)", class = "b", coef = "armND"),
set_prior("normal(0, 5)", class = "b", coef = "condN"),
set_prior("normal(0, 0.01)", class = "sd"),
set_prior("normal(0, 5)", class = "Intercept", dpar = "sigma"),
set_prior("normal(0, 5)", class = "b", dpar = "sigma"),
set_prior("normal(0, 2)", class = "Intercept", dpar = "alpha"),
set_prior("normal(0, 2)", class = "b", dpar = "alpha")
)) prior class coef group resp dpar nlpar lb ub source
normal(0.8, 10) b Intercept <NA> <NA> user
normal(0, 10) b <NA> <NA> user
normal(0, 5) b armND <NA> <NA> user
normal(0, 5) b condN <NA> <NA> user
normal(0, 0.01) sd <NA> <NA> user
normal(0, 5) Intercept sigma <NA> <NA> user
normal(0, 5) b sigma <NA> <NA> user
normal(0, 2) Intercept alpha <NA> <NA> user
normal(0, 2) b alpha <NA> <NA> userif (MODEL) {
m11_transFlat <- brm(m11_transFlat_form,
prior = m11_transFlat_prior,
family = skew_normal(),
init = "0",
data = drum_beats)
m11_transFlat <- add_criterion(m11_transFlat,
"loo",
reloo = TRUE)
save(m11_transFlat,
file = "m11_transFlat.rda")
} else {
load("m11_transFlat.rda")
}Tab. 14. Model summary. The model has not converged.
(m11_transFlat_summary <- summary(m11_transFlat,
priors = TRUE))Warning: There were 8 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: transFlat ~ 0 + Intercept + patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ patient * cond * arm
alpha ~ patient * cond * arm
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
b ~ normal(0, 10)
b_armND ~ normal(0, 5)
b_condN ~ normal(0, 5)
b_Intercept ~ normal(0.8, 10)
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 5)
<lower=0> sd ~ normal(0, 0.01)
Multilevel Hyperparameters:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.01 0.01 0.00 0.02 1.00 3240 2416
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.01 0.01 0.00 0.02 1.00 2715 1855
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.02 0.00 0.02 0.03 1.00 2888 2807
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sigma_Intercept -3.53 0.07 -3.65 -3.40 1.00 2679
alpha_Intercept -4.26 0.74 -5.85 -2.93 1.00 5025
Intercept 0.70 0.02 0.67 0.74 1.00 1817
patient1 0.05 0.02 0.02 0.08 1.00 1548
condN 0.02 0.01 -0.01 0.05 1.00 2225
armND 0.01 0.01 -0.02 0.04 1.00 2364
patient1:condN -0.01 0.00 -0.02 -0.00 1.00 3165
patient1:armND -0.02 0.01 -0.03 -0.01 1.00 2333
condN:armND 0.01 0.01 0.00 0.02 1.00 2554
patient1:condN:armND -0.01 0.01 -0.02 0.00 1.00 2313
sigma_patient1 -0.22 0.10 -0.41 -0.04 1.00 2343
sigma_condN -0.15 0.10 -0.35 0.04 1.00 2417
sigma_armND 0.45 0.09 0.28 0.63 1.00 2253
sigma_patient1:condN -0.16 0.14 -0.44 0.11 1.00 2125
sigma_patient1:armND -0.29 0.13 -0.53 -0.03 1.00 2128
sigma_condN:armND -0.76 0.13 -1.01 -0.52 1.00 1857
sigma_patient1:condN:armND 0.48 0.18 0.12 0.84 1.00 1924
alpha_patient1 0.03 1.05 -2.06 2.02 1.00 5137
alpha_condN 2.37 0.91 0.65 4.18 1.00 4646
alpha_armND -1.38 1.04 -3.48 0.59 1.00 4454
alpha_patient1:condN 0.90 1.24 -1.55 3.34 1.00 4281
alpha_patient1:armND -0.72 1.38 -3.47 1.97 1.00 5135
alpha_condN:armND 0.67 1.24 -1.71 3.12 1.00 4001
alpha_patient1:condN:armND 0.45 1.54 -2.51 3.50 1.00 5173
Tail_ESS
sigma_Intercept 2885
alpha_Intercept 3568
Intercept 2387
patient1 2257
condN 1890
armND 1852
patient1:condN 3400
patient1:armND 2923
condN:armND 3083
patient1:condN:armND 3163
sigma_patient1 3217
sigma_condN 3166
sigma_armND 2540
sigma_patient1:condN 2930
sigma_patient1:armND 2991
sigma_condN:armND 2775
sigma_patient1:condN:armND 2800
alpha_patient1 2958
alpha_condN 3268
alpha_armND 3073
alpha_patient1:condN 2985
alpha_patient1:armND 3109
alpha_condN:armND 2960
alpha_patient1:condN:armND 2845
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(m11_transFlat, ndraws = 100)
Fig. 31. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
beautify_my_plot(plot(conditional_effects(m11_transFlat,
dpar = "mu"),
plot = FALSE))
conditions <- make_conditions(drum_beats,
vars = "arm")
beautify_my_plot(plot(conditional_effects(m11_transFlat,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11_transFlat,
dpar = "sigma"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11_transFlat,
dpar = "sigma",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11_transFlat,
dpar = "alpha"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11_transFlat,
dpar = "alpha",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
Fig. 32. Conditional plots. Had to take out all the fancy figure polishing, subplot lettering etc. Subplots are supposed to be numbered from A to … starting in the top left corner, proceeding row-wise.
2.2.2 Conclusions
wrap_plots(plot(conditional_effects(m11_transFlat,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
plot = FALSE)[[1]] +
scale_color_colorblind() +
theme(legend.position = "none"),
plot(conditional_effects(m11_transFlat,
dpar = "sigma",
effects = "patient:cond",
conditions = conditions),
plot = FALSE)[[1]] +
scale_color_colorblind() +
theme(legend.position = "none"),
plot(conditional_effects(m11_transFlat,
dpar = "alpha",
effects = "patient:cond",
conditions = conditions),
plot = FALSE)[[1]]) +
scale_color_colorblind() +
plot_annotation(tag_levels = "A")
Fig. 30. Three-way interaction plots of the three modeled parameters \(\mu\), \(\sigma\), and \(\alpha\) as estimated by the full model. The x axis shows group membership (0 = healthy subjects, 1 = patients). The instruction conditions are color-coded (black = controlled strokes, brown = normal strokes), and the facets show the side (1 = dominant, 2 = non-dominant).
In transFlat there is a small but most likely positive
difference between controls and patients, but all the other main effects
are not clearly different from zero
(Fig. 30A).
Variation shows a curiously extreme estimated value for healthy
subjects playing controlled strokes on the non-dominant side
(Fig. 30B) which seems to be
mostly driven by a relatively large estimated value for the non-dominant
side (sigma_armND in Tab.
14).
Fig. 30C shows that the left
tail of transFlat’s distribution gets heavier in the
controlled condition, suggesting a tendency in players to produce a
larger amount of low-valued transFlat values when
performing unfamiliar strokes.
Since there is no clear difference in the estimated location
parameter \(\mu\) for
transFlat values between conditions, this descriptor is
probably not a suitable cue to distinguish between conditions in a
listening test.
2.3 Univariate Models –
transCrest
2.3.1 The Full Model: Modeling the Interactions of Group, Instruction, and Side on Shape and Spread
(m11_transCrest_form <-bf(transCrest ~ 0 + Intercept + patient * cond * arm +
(1 | subj) +
(1 | cond) +
(1 | arm),
sigma ~ patient * cond * arm,
alpha ~ patient * cond * arm))transCrest ~ 0 + Intercept + patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ patient * cond * arm
alpha ~ patient * cond * arm(m11_transCrest_prior <- c(set_prior("normal(2.5, 10)", class = "b", coef = "Intercept"),
set_prior("normal(0, 10)", class = "b"),
set_prior("normal(0, 5)", class = "b", coef = "armND"),
set_prior("normal(0, 5)", class = "b", coef = "condN"),
set_prior("normal(0, 0.01)", class = "sd"),
set_prior("normal(0, 5)", class = "Intercept", dpar = "sigma"),
set_prior("normal(0, 5)", class = "b", dpar = "sigma"),
set_prior("normal(0, 2)", class = "Intercept", dpar = "alpha"),
set_prior("normal(0, 2)", class = "b", dpar = "alpha")
)) prior class coef group resp dpar nlpar lb ub source
normal(2.5, 10) b Intercept <NA> <NA> user
normal(0, 10) b <NA> <NA> user
normal(0, 5) b armND <NA> <NA> user
normal(0, 5) b condN <NA> <NA> user
normal(0, 0.01) sd <NA> <NA> user
normal(0, 5) Intercept sigma <NA> <NA> user
normal(0, 5) b sigma <NA> <NA> user
normal(0, 2) Intercept alpha <NA> <NA> user
normal(0, 2) b alpha <NA> <NA> userif (MODEL) {
m11_transCrest <- brm(m11_transCrest_form,
prior = m11_transCrest_prior,
family = skew_normal(),
init = "0",
control = list(max_treedepth = 15),
data = drum_beats)
m11_transCrest <- add_criterion(m11_transCrest,
"loo",
reloo = TRUE)
save(m11_transCrest,
file = "m11_transCrest.rda")
} else {
load("m11_transCrest.rda")
}Tab. 15. Model summary.
(m11_transCrest_summary <- summary(m11_transCrest,
priors = TRUE)) Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: transCrest ~ 0 + Intercept + patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm)
sigma ~ patient * cond * arm
alpha ~ patient * cond * arm
Data: drum_beats (Number of observations: 1102)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Priors:
b ~ normal(0, 10)
b_armND ~ normal(0, 5)
b_condN ~ normal(0, 5)
b_Intercept ~ normal(2.5, 10)
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 5)
<lower=0> sd ~ normal(0, 0.01)
Multilevel Hyperparameters:
~arm (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.01 0.01 0.00 0.02 1.00 3013 1546
~cond (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.01 0.01 0.00 0.02 1.00 3522 1983
~subj (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.09 0.01 0.08 0.10 1.00 4434 2804
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sigma_Intercept -0.96 0.06 -1.08 -0.83 1.00 2792
alpha_Intercept 5.57 0.91 3.94 7.53 1.00 3944
Intercept 3.03 0.06 2.92 3.15 1.00 2042
patient1 -0.51 0.08 -0.66 -0.36 1.00 1965
condN -0.15 0.04 -0.23 -0.08 1.00 2481
armND -0.18 0.05 -0.27 -0.09 1.00 2511
patient1:condN 0.08 0.05 -0.01 0.16 1.00 2428
patient1:armND 0.30 0.06 0.18 0.41 1.00 2313
condN:armND 0.10 0.06 -0.02 0.22 1.00 2326
patient1:condN:armND -0.18 0.07 -0.33 -0.04 1.00 2269
sigma_patient1 -0.43 0.10 -0.62 -0.24 1.00 2580
sigma_condN -0.63 0.09 -0.80 -0.45 1.00 2853
sigma_armND 0.05 0.09 -0.12 0.22 1.00 2428
sigma_patient1:condN 0.49 0.14 0.21 0.75 1.00 2631
sigma_patient1:armND 0.42 0.14 0.15 0.69 1.00 2213
sigma_condN:armND 0.63 0.12 0.39 0.87 1.00 2378
sigma_patient1:condN:armND -0.91 0.18 -1.26 -0.54 1.00 2342
alpha_patient1 -0.06 1.17 -2.32 2.24 1.00 3764
alpha_condN -3.19 0.98 -5.15 -1.29 1.00 3858
alpha_armND 1.04 1.26 -1.43 3.53 1.00 3746
alpha_patient1:condN 0.54 1.19 -1.87 2.88 1.00 3727
alpha_patient1:armND 3.86 1.51 0.98 6.80 1.00 4501
alpha_condN:armND 0.22 1.34 -2.37 2.92 1.00 3836
alpha_patient1:condN:armND 1.51 1.71 -1.75 4.86 1.00 5090
Tail_ESS
sigma_Intercept 3223
alpha_Intercept 3185
Intercept 2562
patient1 2268
condN 2898
armND 2964
patient1:condN 2827
patient1:armND 2950
condN:armND 2956
patient1:condN:armND 2671
sigma_patient1 3033
sigma_condN 3190
sigma_armND 2700
sigma_patient1:condN 2689
sigma_patient1:armND 2574
sigma_condN:armND 3004
sigma_patient1:condN:armND 2649
alpha_patient1 2774
alpha_condN 3036
alpha_armND 2685
alpha_patient1:condN 2809
alpha_patient1:armND 2981
alpha_condN:armND 2799
alpha_patient1:condN:armND 3189
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(m11_transCrest, ndraws = 100)
Fig. 34. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.
beautify_my_plot(plot(conditional_effects(m11_transCrest,
dpar = "mu"),
plot = FALSE))
conditions <- make_conditions(drum_beats,
vars = "arm")
beautify_my_plot(plot(conditional_effects(m11_transCrest,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11_transCrest,
dpar = "sigma"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11_transCrest,
dpar = "sigma",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11_transCrest,
dpar = "alpha"),
plot = FALSE))
beautify_my_plot(plot(conditional_effects(m11_transCrest,
dpar = "alpha",
effects = "patient:cond",
conditions = conditions),
plot = FALSE))
Fig. 35. Conditional plots. Had to take out all the fancy figure polishing, subplot lettering etc. Subplots are supposed to be numbered from A to … starting in the top left corner, proceeding row-wise.
2.3.2 Conclusions
wrap_plots(plot(conditional_effects(m11_transCrest,
dpar = "mu",
effects = "patient:cond",
conditions = conditions),
plot = FALSE)[[1]] +
scale_color_colorblind() +
theme(legend.position = "none"),
plot(conditional_effects(m11_transCrest,
dpar = "sigma",
effects = "patient:cond",
conditions = conditions),
plot = FALSE)[[1]] +
scale_color_colorblind() +
theme(legend.position = "none"),
plot(conditional_effects(m11_transCrest,
dpar = "alpha",
effects = "patient:cond",
conditions = conditions),
plot = FALSE)[[1]]) +
scale_color_colorblind() +
plot_annotation(tag_levels = "A")
Fig. 30. Three-way interaction plots of the three modeled parameters \(\mu\), \(\sigma\), and \(\alpha\) as estimated by the best-fitting model. The x axis shows group membership (0 = healthy subjects, 1 = patients). The instruction conditions are color-coded (black = controlled strokes, brown = normal strokes), and the facets show the side (1 = dominant, 2 = non-dominant).
The location parameter \(\mu\) of
the transCrest data has been estimated to be clearly
different between groups, conditions, and sides such that there is even
a three-way interaction between the main effects
(Fig. 30A): patients have low
\(\mu\) on the dominant side and higher
\(\mu\) on the non-dominant side, while
the opposite holds for healthy controls.
The estimation of spread displays some weird patterns, with \(\sigam\) generally being larger in the in the controlled instruction condition, but with one strangely small estimate in patients on the dominant side. The opposite holds for the normal instruction condition, with a saliently large value in healthy subjects when playing normal strokes with the non-dominant arm.
The distributio of the data has a heavy tail on the right side seems to follow a pattern of “normal” < “controlled”, and “dominant” < “non-dominat”. Only for patients on the non-dominant side skewness seems to be about the same, irrespective of instruction.
With the three-way interaction in the location parameter this
descriptor seems a promising candidate for a listening test because if
we aim for a single descriptor to be sufficient for listeners to
reliably distinguish between all three factors, only a descriptor
satisfying this interaction can be used to delineate the different
strokes. But since the original aim–as I understand it–is to primarily
distinguish between controlled and normal strokes, every descriptor with
a clear difference between these two factor levels would be a good
candidate from a statistical point of view. In the three descriptors
investigated here, transSC most clearly distinguishes
(again: from a stats point of view) between normal and controlled
instruction conditions.
3 System Setup
Data wrangling and analyses were carried out with the statistical
package R [R version 4.4.2 (2024-10-31); R Core Team (2020)]. Bayesian modeling was done
with the package brms (Bürkner 2017,
2018) which uses the probabilistic language Stan as back end
(Carpenter et al. 2017). Plots were done
with the packages bayesplot (Gabry
and Mahr 2020) and ggplot2 (Wickham 2016).